Data fitting

Curve fitting is the process of constructing a curve, or mathematical function, that has the best fit to a series of eperimental data points.

In this exercise we will calculate the sum of the "medipix.edf" image along the vertical and horizontal axis and fit a gaussian function, either using scipy or the module from PyMca.

Initialization of the scientific scripting environment


In [1]:
%pylab inline


Populating the interactive namespace from numpy and matplotlib

In [2]:
import fabio
img = fabio.open("medipix.edf").data
imshow(np.log(img))


-c:3: RuntimeWarning: divide by zero encountered in log
Out[2]:
<matplotlib.image.AxesImage at 0x7ff5f78fb2d0>

In [3]:
vert = img.sum(axis=1)
plot(vert)
title("Sum along the horizontal direction")


Out[3]:
<matplotlib.text.Text at 0x7ff5f78fbc10>

Fitting function definition

The fitting function is a function of the data-point (xdata) and of a sert of parameters.

For a gaussian function defined as: gaussian


In [4]:
def gaussian(x, height, center, sigma, offset):
    return height*exp(-(x-center)**2/(2*sigma**2)) + offset

An essential part is the initial set of parameters p0:


In [5]:
p0 = [70000, 140, 10,10000]
xdata = numpy.arange(len(img),dtype="float")
ydata = vert.astype("float")

In [6]:
from scipy.optimize import curve_fit
p,cov = curve_fit(gaussian, xdata, ydata, p0)
print(p)


[  6.31031426e+04   1.35955545e+02   2.99885986e+00   1.19186952e+04]

In [7]:
#Display the result of the fit
plot(gaussian(xdata, *p0), label="Initial guess")
plot(gaussian(xdata, *p), label="Fitted curve")
plot(vert,label="Experimental data")
legend()
title("Result of the fit")


Out[7]:
<matplotlib.text.Text at 0x7ff5f3c89510>

Using PyMca library


In [8]:
try:
    from PyMca5.PyMca import Gefit
except ImportError:
    from PyMca import Gefit

PyMca's Levenberg-Marquardt fitting module needs the function to be layed out the other way around:


In [9]:
def gaussian_ge(param, xdata):
    return gaussian(xdata, *param)

p, chi2, err = Gefit.LeastSquaresFit(gaussian_ge, p0, xdata=xdata, ydata=ydata)
print(p)


[63103.22036555503, 135.95551675017475, 2.998849240794874, 11918.81798431471]

In [10]:
#Display the result of the fit
plot(gaussian(xdata, *p0), label="Initial guess")
plot(gaussian(xdata, *p), label="Fitted curve")
plot(vert,label="Experimental data")
legend()
title("Result of the fit")


Out[10]:
<matplotlib.text.Text at 0x7ff5f3bdcb10>

Conclusion

Curve fitting goes alway via the following steps:

  • get the data_point x_data and y_data
  • define the fitting function as y_data = function (x_data, param)
  • chose an initial guess for the set of parameters p0
  • run the optimizer
  • check the result.

Note that some optimizer don't fit the function but instead minimize the error_function: error(param) = y_data - function(x_data) It is often necessary to look at the documentation of the fitting function


In [ ]: